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Abstract. We examine the cosmological constraining power of future large-scale weak lensing surveys on 
the model of the ESA planned mission Euclid, with particular reference to primordial non-Gaussianity. Our 
analysis considers several different estimators of the projected matter power spectrum, based on both shear and 
flexion. We review the covariance and Fisher matrix for cosmic shear and evaluate those for cosmic flexion 
and for the cross-correlation between the two. The bounds provided by cosmic shear alone are looser than 
previously estimated, mainly due to the reduced sky coverage and background number density of sources for 
the latest Euclid specifications. New constraints for the local bispectrum shape, marginalized over erg, are at 
the level of A/nl ~ 100, with the precise value depending on the exact multipole range that is considered in the 
analysis. We consider three additional bispectrum shapes, for which the cosmic shear constraints range from 
A/nl ~ 340 (equilateral shape) up to A/nl ~ 500 (orthogonal shape). Also, constraints on the level of non- 
Gaussianity and on the amplitude of the matter power spectrum cr 8 are almost perfectly anti-correlated, except 
for the orthogonal bispectrum shape for which they are correlated. The competitiveness of cosmic flexion 
constraints against cosmic shear ones depends by and large on the galaxy intrinsic flexion noise, that is still 
virtually unconstrained. Adopting the very high value that has been occasionally used in the literature results in 
the flexion contribution being basically negligible with respect to the shear one, and for realistic configurations 
the former does not improve significantly the constraining power of the latter. Since the shear shot noise is 
white, while the flexion one decreases with decreasing scale, by considering high enough multipoles the two 
contributions have to become comparable. Extending the analysis up to ^ max = 20, 000 cosmic flexion, while 
being still subdominant, improves the shear constraints by ~ 10% when added. However on such small scales 
the highly non-linear clustering of matter, the impact of baryonic physics, and the non-Gaussian part of the 
covariance matrix make any error estimation uncertain. By considering lower, and possibly more realistic, 
values of the flexion intrinsic shape noise results in flexion constraining power being a factor of ~ 2 better than 
that of shear, and the bounds on cr 8 and /nl being improved by a factor of ~ 3 upon their combination. 
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1 Introduction 

In the forthcoming decade measurements of the fundamental parameters of cosmology will undergo a signif- 
icant enhancement in precision, especially thanks to a number of space-based missions. The satellite Planck 
[1] has been taking data of the Cosmic Microwave Background (CMB henceforth) temperature and polariza- 
tion fluctuations for about two years now, and their cosmological interpretation should soon be available. The 
X-ray observatory eROSITA [2] is scheduled for launch in 2013, and will provide an all-sky catalog of galaxy 
clusters with unprecedent resolution and sensitivity for a survey of this size. Such a vast X-ray cluster sample 
is ideally suited for cosmology, as the uncertainties in the relations between cluster mass and observables can 
be treated as 'nuisance' parameters and marginalized over [3]. In the more distant future, the Euclid 1 satellite 
[4] is currently scheduled for launch in 2019, providing numerous opportunities for cosmological investiga- 
tion with optical/NIR data: spectroscopy of a large number of galaxies over the almost entire extragalactic 
sky will allow measurements of the Baryon Acoustic Oscillation (BAO); spectroscopically, photometrically, 
and weak lensing-selected clusters of galaxies will provide a massive sample ideal for cosmological parameter 
estimation; the weak lensing map of almost half the sky will chart the matter power spectrum over unprecedent 
ranges of scales and redshifts, thus providing fundamental clues on the growth of matter perturbations. Planned 
by NASA, the Wide-Field InfraRed Survey Telescope (WFIRST, [5]) is supposed to have mission objectives 
similar to those of Euclid, although with more focus on supernova cosmology. 

1 http://www.euclid-ec.org 
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Given the enormous quantity of data that each one of these projects is expected to return, as well as 
the vast amount of resources that are invested, it is utterly important to explore different ways in which these 
products can be combined, and establish which data or data combinations are most promising for a certain 
cosmological application. This motivates us to review in this paper the constraining power of one of the most 
immediate science products for the Euclid mission, that is the power spectrum of galaxy shapes. As galaxy 
shapes are affected by the light deflection operated by the Large-Scale Structure (LSS), their power spectrum 
is basically a projected (and weighted) version of the three-dimensional non-linear matter power spectrum, 
thus containing a wealth of cosmological information. The moment of the galaxy light distribution that is most 
commonly used is the ellipticity, since it provides a measure of the shear field [6]. However, higher order 
lensing fields, affecting different light moments, exist and have been measured for individual galaxy clusters 
[7, 8]. Some of these fields, such as the flexion fields [9], will be measurable on cosmological scales by future 
weak lensing surveys such as Euclid and WFIRST. In this work we assess the level at which cosmic shear and 
flexion, separately and in combination, can be used in order to constrain Primordial Non-Gaussianity (PNG), 
adopting Euclid as our fiducial weak lensing survey. 

According to the standard picture for structure formation, the primordial density fluctuations that seeded 
today's LSS were produced during an early accelerated expansion phase of the Universe, dubbed inflation 
[10]. The simplest inflationary model predicts almost Gaussian density fluctuations, while numerous more 
sophisticated models forecast arbitrary deviations from Gaussianity. If occurring, PNG would have a significant 
effect both on CMB and cosmic structures. As a matter of fact, the issue of constraining deviations from 
primordial Gaussianity by using structure formation and evolution has recently attracted much attention in 
the literature, with efforts directed towards the abundance of dark matter halos [11-16], halo biasing [17- 
22]), galaxy bispectrum [23, 24], mass density distribution [25] and topology [26, 27], integrated Sachs-Wolfe 
effect [28-30], Lyor flux from low-density intergalactic medium [31], 21-cm fluctuations [32, 33], reionization 
[34], and weak lensing peak counts [35, 36]. PNG is one of the very few handle that we have on the pre- 
recombination universe, and it is hence important to study it in detail and with different probes. 

Cosmic shear from future weak lensing surveys has been investigated as a tool for probing PNG in [37] 
(see also [38]). Here we review those constraints with a more complete statistical analysis and by using the 
latest Euclid specifications for the galaxy number density, the intrinsic ellipticity noise, the sky coverage, and 
the source redshift distribution. In addition, we consider cosmic flexion as an alternative channel of cosmo- 
logical information. The intrinsic shape noise for flexion has not been reliably measured yet. Values adopted 
in the literature are usually very high, hence they are expected to result in flexion giving only a subdominant 
contribution to the overall Euclid constraining power. However, given the large uncertainty on the flexion noise 
and the fact this noise is scale-dependent, it is important to investigate if the flexion contribution is indeed neg- 
ligible and if there are configurations for which it can become relevant. We dedicate the Appendix A to the 
calculations of covariances and Fisher matrices for cosmic shear, cosmic flexion, and their cross-correlation, 
the latter two of which, to the best of our knowledge, have not been reported elsewhere. The main body of this 
paper is organized as follows. In Section 2 we review how PNG is parametrized and in Section 3 we summarize 
the effect thereof on halo mass function, bias, and internal structure. In Section 4 we describe the halo model, 
a physically motivated framework that we adopted in order to estimate the fully non-linear matter and lensing 
power spectra in non-Gaussian cosmologies. In Section 5 we show our results concerning joint constraints 
on the level of PNG and the amplitude of the matter power spectrum from cosmic shear and cosmic flexion. 
Section 6 is dedicated to a discussion of the results and a summary of the conclusions. 

Throughout this work we adopted as a fiducial cosmological model the one suggested by the latest anal- 
ysis of the WMAP data [39], bearing Q m , = 0.272, O A , = 1 - O m>0 , Q b ,o = 0.046, H = /il00 km s" 1 Mpc 1 
with h = 0.704, and a matter power spectrum normalized by setting cr g = 0.809. 

2 Non- Gaussian cosmologies 
2.1 General 

As mentioned in the Introduction, extensions of the most standard model of inflation [10, 40, 41] can produce 
substantial deviations from a Gaussian distribution of primordial density and potential fluctuations (see [42, 43] 
for recent reviews). The amount and shape of these deviations depend critically on the kind of non-standard 
inflationary model that one has in mind, as will be detailed later on. 
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A particularly convenient (although not unique) way to describe generic deviations from Gaussian initial 
conditions consists in writing the gauge-invariant Bardeeen's potential <1> as the sum of a Gaussian random field 
and a quadratic correction [12, 44-^-6], according to 



The parameter /nl in Eq. (2.1) determines the amplitude of non-Gaussian deviations, and it is in general depen- 
dent on the scale. The symbol * denotes convolution between functions, and reduces to standard multiplication 
upon constancy of /nl- In the following we adopted the LSS convention (as opposed to the CMB convention, 
see [15, 28, 29, 47]) for the definition of the fundamental parameter /nl- According to this, the primordial 
value of <t> has to be linearly extrapolated at z = 0, and as a consequence the constraints given on /nl by the 
CMB have to be raised by ~ 30 per cent to comply with this paper's convention (see also [20] for a concise 
explanation). 

If /nl # the potential i> is a random field with a non-Gaussian probability distribution. Therefore, the 
field itself cannot be described by the linear power spectrum Po.l(^) = Ak"~ 4 alone, rather higher order statis- 
tics are needed. In many circumstances the dominant higher order statistic is the bispectmm Bq,{k\,k2, k 3 ). 
The bispectrum is the Fourier transform of the three-point correlation function (O^OO^)®^)) and it can 
hence be implicitly defined as 



where 6d is the Dirac delta distribution and the angular brackets indicate an ensemble average. 

Since different inflationary models predict different non-Gaussian shapes, that is different ways in which 
the potential bispectrum depends on its three arguments, understanding the PNG shape is of fundamental 
importance in order to pinpoint the physics of the early Universe. In the standard inflationary scenario the 
early accelerated expansion is accounted for by a slowly-rolling scalar field, the inflaton. In the present work 
we considered four different shapes of the potential bispectrum, arising from different modifications of this 
standard scenario. They are all genetically described in the following. For further details we refer the reader to 
[21]. 

2.2 Shapes 

2.2.1 Local shape 

The standard single-field inflationary scenario generates negligibly small deviations from Gaussianity. These 
deviations are commonly said to be of the local shape, and the related bispectrum of the Bardeen's potential 
is maximized for squeezed configurations, where one of the three wavevectors has a much smaller magnitude 
than the other two. In this case the parameter /nl must be a constant, and it is expected to be of the same order 
of the slow-roll parameters [48], that are in turn of order unity. 

However non-Gaussianities of the local shape can also be generated in the case in which an additional 
light scalar field, different from the inflaton, contributes to the observed curvature perturbations [49]. This 
happens, for instance, in curvaton models [50, 51] or in multi-field models [52, 53]. In this case the parameter 
/nl is allowed to be substantially different from zero. The best constraints on the level of non-Gaussianity for 
a local bispectrum shape come from the WMAP-7 data [39], and at 95% Confidence Level (CL) amount to 
-10 < / NL < 74. 

2.2.2 Equilateral shape 

In some inflationary models the kinetic term of the inflaton Lagrangian is not standard, containing higher- 
order derivatives of the field itself. One significant example of this is the DBI model ([54, 55], see also 
[56-58]). In this case the primordial bispectrum is maximized for configurations where the three wavevectors 
have approximately the same amplitude. A template for this equilateral bispectrum can be found in [59]. 

Unlike the local shape, for the equilateral case there is no theoretical prescription against a running of the 
/nl parameter with the scale, and in fact many authors introduce such a running in order to enhance deviations 
from Gaussianity at scales smaller than those probed by CMB studies. In this case however we opted for not 
introducing such a running, in order to directly compare pure bispectrum shapes. Constraints on /nl in this 
case have also been given by the WMAP team [39], and correspond to -214 < /nl < 266 at 95% CL. 




(2.1) 



(®(ki)<S>(k 2 )®(h)) = (2^) j 5 D (A: 1 + k 2 + k 3 )B^k u k 2 ,k 3 ) , 



(2.2) 
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2.2.3 Enfolded shape 

For deviations from Gaussianity evaluated in the regular Bunch-Davies vacuum state, the primordial potential 
bispectrum is of local or equilateral shape, depending on whether or not higher-order derivatives play a signifi- 
cant role in the evolution of the inflaton field. If the Bunch-Davies vacuum hypothesis is dropped, the resulting 
bispectrum is maximal for squashed configurations [60, 61]. Authors in [62] found a template that describes 
very well the properties of this enfolded-shape bispectrum. 

Constraints on the level of PNG for the enfolded-shaped bispectrum have not been given in [39] . However, 
relatively tight bounds can be estimated from the correlation functions of LSS tracers. Specifically, authors in 
[63] give -12 < / NL < 358 (CMB convention) at the 2cr CL. 

2.2.4 Orthogonal shape 

A shape of the bispectrum can be constructed that is nearly orthogonal (with respect to a suitably defined scalar 
product) to both the local and equilateral forms [64]. This kind of bispectrum is peaked both on equilateral and 
squashed configurations, with opposite signs. Constraints on the level of non-Gaussianity compatible with the 
CMB in the orthogonal scenario were also given by the WMAP team [39], corresponding to -410 < /nl < 6 
at 95% CL. 

Similarly to the equilateral case, for enfolded and orthogonal shapes it would be allowed to introduce a 
running of the parameter /nl with scale, which we decided not to include. However it should be stressed that, 
differently from the equilateral shape, in the two latter cases there is no first principle that can guide one in the 
choice of a particular kind of running, and until now no work has addressed the problem of a running for these 
shapes [65, 66]. 



3 Impact of PNG on structure formation 

PNG produces modifications in the statistics of density peaks, resulting in differences in the mass function, 
bias, and internal structure of dark matter halos. Since all these factors enter in the modeling of the LSS, we 
summarize in the following how these modifications can be parametrized. 

3.1 Mass function 

We modeled the non-Gaussian contribution to the number counts of dark matter halos according to the pre- 
scription presented in [14]. These authors gave a useful expression for the Press & Schechter [67] mass func- 
tion following from non-Gaussian initial conditions, «ps(M, z). Since the corresponding mass function com- 
puted with Gaussian initial conditions, nj?(M, z), is well known, we can define a correction factor fi(M, z) = 
«ps(M, z)/«pg'(M, z). Under the assumption that the effect of PNG on the mass function is independent of 
the prescription adopted to describe the mass function itself, it follows that the non-Gaussian mass function 
computed according to an arbitrary prescription, «(M, z), can be related to its Gaussian counterpart through 

n(M, z) = «(M, z) n (G) (M, z) . (3.1) 

In order to evaluate «ps(M, z), and hence ^(M, z), authors in [14] performed an Edgeworth expansion 
[68] of the probability distribution for the smoothed density fluctuations field, truncating it at the linear term 
in cr M . One function that turns out to be needed in the subsequent formula is the reduced skewness Si{M) = 
/nl f*?,(M)/o- 4 M of the non-Gaussian distribution, where the skewness fi3(M) can be computed as 

{M)= r d 3 fcid 3 fc 2 d 3 fc 3 MRikdMRi^MRihXQihmkimm . (3.2) 

Jr* (2?r) y 

The function Mr(Ic) relates the density fluctuations smoothed on some scale R (corresponding to the mass M) 
to the respective peculiar potential, 

M R (k) = l T( f ecl W R (k) , (3.3) 
3 ff2 Qm0 
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where T(k) is the matter transfer function and WrQi) is the Fourier transform of the top-hat window function. 
In this work we adopted the Bardeen et al. [69] matter transfer function, with the shape factor correction 
presented in [70]. This reproduces fairly well the more sophisticated recipe of Eisenstein & Hu [71] except for 
the presence of the BAO, that anyway is not of interest here. 

For the reference Gaussian mass function we adopted the Sheth & Tormen prescription [72] (see [73- 
75] for alternative prescriptions). It is worth noting that other approaches for computing the non-Gaussian 
correction to the halo number counts exist (see for instance [11]). These alternative approaches give results that 
are in broad agreement with those of [14], and have been collectively shown to provide a good representation 
of the mass function measured in cosmological simulations [15, 76]. 

3.2 Halo bias 

In this Subsection we describe how the bias of dark matter halos gets modified by PNG. The fact that deviations 
from primordial Gaussianity induce a scale-dependence on the large-scale halo bias has been first noticed by 
[17]. Later, authors in [18] used the peak-background split formalism in order to provide a semi-analytic 
expression for the non-Gaussian bias, hence allowing to place and forecast stringent constraints on PNG by 
using the correlation functions of LSS tracers [28, 29, 63, 77, 78]. 

The non-Gaussian halo bias can be written in a relatively straightforward way in terms of its Gaussian 
counterpart as [79] 



where the function /?«(£) encapsulates all the scale dependence of the non-Gaussian correction to the bias. It 
can be written as 



with a 2 = k 2 + ( 2 + 2k(fi. In the simple case of local bispectrum shape it can be shown that the function 
should scale as oc k~ 2 at large scales, so that a substantial boost (for positive /nl) in the halo bias is expected at 
those scales. For the Gaussian halo bias we adopted the Sheth, Mo & Tormen [80] prescription. 

We emphasize that the semi-analytic prescription by [18] has also been extensively tested against numer- 
ical simulations, giving an overall good agreement [15, 81]. 

3.3 Halo internal structure 

Since PNG affects the timing of structure formation, it is expected to have some sort of impact on the internal 
structure of dark matter halos. However, to date there are only a handful of works exploring this matter. 
After the pioneering study by [82], more recently this issue has been addressed in passing by [83] through 
numerical simulations, and by [84] via semi-analytic considerations. These authors assume the average dark 
matter halo to be well represented by a Navarro, Frenk & White [85] (NFW henceforth) density profile also in 
non-Gaussian cosmologies. However the concentration of the profile turns out to be enhanced in models with 
positive skewness and depressed in models with negative skewness. The results of the semi-analytic model of 
[84] fairly agree with the simulations by [83], showing that the halo concentration is increased (decreased) by 
~ 4 — 10% for a local bispectrum shape with /nl = +100 (/nl = -100, CMB convention). The effect is mildly 
dependent on mass and redshift, with more massive objects placed at higher redshift being the most affected by 
PNG. The semi-analytic model by [84] is however quite cumbersome. A simplified model has been proposed 
by [86]. The latter reproduces well the results presented in [84] at low redshift and on galaxy group-scales, but 
it somewhat underestimates them for large masses/high redshifts. 

In this work we assumed average dark matter halos to be well represented by NFW density profiles, 
however we adopted the same concentration prescription for all cosmological models, which matches the 
results of ACDM cosmological simulations (see [37] for more details). The reason for this choice is twofold: 
;) the previous studies are limited to the local bispectrum shape, while in this work we explored a variety of 
non-Gaussian shapes. In principle the model by [84] is of general application, however its agreement with 
cosmological simulations has not been tested so generally; if) in general we limited our analysis to relatively 
large scales, where the impact of the very inner structure of dark matter halos can be considered irrelevant. 




(3.4) 




(3.5) 
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In the few cases in which it is not (that will be clearly indicated), one should keep in mind that the impact of 
baryonic physics and the non-Gaussian parts of the covariances (see Appendix A for details) constitute more 
of an uncertainty in the statistical analysis. We return on this point further below. 



4 Modeling the large-scale structure 
4.1 Halo model 

We represented the non-linear matter power spectrum by making use of the halo model [87, 88]. This is a 
physical framework that allows the description of the correlation of dark matter particles as well as that of 
different tracers of the LSS such as galaxies ad galaxy clusters. It is based on the fact that the power spectrum 
of particles (either dark matter particles or tracers) is given by the sum of two contributions: particle pairs 
residing in the same structure, and particle pairs residing in two different structures. Accordingly, the dark 
matter power spectrum can be written as the sum of two contributions, 



P m (k,z) = P m ,i(k,z) + P m , 2 (k,z). (4.1) 

The first term is named 1-halo term, and dominates on very small scales, while the second is the 2-halo term, 
and it is dominant at large scales. It is easy to understand that the 2-halo term should incorporate the bias of 
dark matter halos, since it represents the clustering of particle pairs residing in separated structures [89]. 
The two contributions to the power spectrum of dark matter particles can be written as 



dM n(M,z) 

o 

and 



p(M,z,k)f 



(4.2) 



P m .i(k,z) 



Jr»+oo 




dM n(M,z)b{M,z,k)^ 



P m ,L(k,z). (4.3) 



Jo Pm 

In the previous pair of equations, p m represents the comoving matter density, Pm^ik, z) is the linear matter 
power spectrum, n(M, z) is the mass function, and b(M, z, k) is the halo bias, where we inserted an additional 
scale dependence in order to account for the effect of primordial non-Gaussianity (see Section 3.2). The 
function p(M, z, k) represents the Fourier transform of the average density profile of dark matter halos, 



f 

Jo 



p(M,z,k)= I r z drp(M,z,r)— — - . (4.4) 

kr 



Since for p(M, z, k) we adopted a NFW shape, its Fourier transform can be computed analytically as described 
in [90, 91]. Because the integral in the previous equation is evaluated up to the virial radius R v of the structure, 
it is easy to see that in the large scale limit p(M, z, k) converges to the virial mass of the halo. 

There are two assumptions underlying the halo model representation of the dark matter power spectrum. 
The first one is that all the matter in the Universe is contained within halos of some mass, that is 

p + OO 

I MdMn(M,z) = p m . (4.5) 
Jo 

As can be verified, this assumption is fulfilled by the Press & Schechter [67] and Sheth & Tormen [72] mass 
functions, of which we are considering the latter. The second assumption is that the fully non-linear power 
spectrum should converge to the linear one at very large scales. In other words, this implies the non-trivial 
constraint 

Jr»+oo 
MdM n(M, z)b(M, z, k) = p m (4.6) 
o 

for very large scales. This constraint has to be explicitly enforced, since the standard halo bias formulations 
do not satisfy it automatically. For this and other details on the practical implementation of the halo model we 
refer the interested reader to [37]. 
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As a last note, in [83] the authors pointed out that a better agreement of the non-linear power spectrum 
estimated by the halo model with numerical simulations may be obtained by considering the so-called halo 
exclusion term. This consists in an integral to be subtracted by the 2-halo contribution, taking into account 
that dark matter halos are considered as hard spheres that cannot overlap with one another. We ignored this 
possible improvement, as what matters in our statistical analysis is only the relative change induced by PNG 
on the non-linear matter power spectrum. 

4.2 Cosmological weak lensing 

An effective way to probe the large-scale matter distribution in the Universe is offered by the gravitational 
deflection of light. Specifically, the weak distortion of images of background galaxies as their light passes 
through the LSS allows to map the total matter distribution up to the source redshift, projected along the line 
of sight. Therefore, the power spectrum of lensing-induced distortions can be thought as a projection of the 
matter power spectrum described in the previous Subsection, weighted with a function that depends on the 
redshift distribution of background sources. 

The quantity relevant for describing the correlation function of lensing-induced distortions is the conver- 
gence power spectrum, that can be related to the fully non-linear matter power spectrum through [6, 92] 



where x - X(z) is the comoving distance out to redshift z, the scale factor a(x) is normalized such that a(0) = 
1, and /k(x) i s me comoving angular diameter distance, which in general depends on the curvature of the 
Universe. In this work we considered flat universes only, so that /k(x) - however we keep the formalism as 
general as possible. As mentioned, the weight function W(x) is related to source redshift distribution g(z) by 



Both the integrals above formally extend up to the comoving horizon distance Xh, however the source redshift 
distribution declines to zero well before that, so that the integral can be effectively thought as extending out to 
infinity. 

It should be kept in mind that the previous equation for the convergence power spectrum is not generally 
applicable, rather it makes use of the Limber approximation for flat-sky quantities. This approximation is valid 
only as long as we consider sufficiently small patches of the sky. Specifically, it has been shown by [24] that 
the accuracy of the Limber approximation for the convergence power spectrum is at the level of 1 % for £ > 10, 
corresponding to angular scales smaller than 2 x 10 3 arcmin. 

The convergence power spectrum cannot be measured directly, but it can be estimated by measuring the 
correlation functions of other (observable) lensing fields. The field that has been mainly used in the past for 
this purpose is the shear y, because it affects only the ellipticity of background sources and it is hence relatively 
easy to measure. The two flexion fields, which affect higher-order moments of the galaxy light distribution have 
been measured recently in galaxy cluster fields, and will certainly be possible to measure them on cosmological 
scales with Euclid. In this work we ignored the second flexion because it is trickier to measure, and focused 
instead only on the first flexion ip, which induces an arc-like distortion on background sources [93]. In the 
Appendix A we summarize how unbiased estimators of the convergence power spectrum can be constructed 
which are based on both shear (cosmic shear) and flexion (cosmic flexion) measurements, as well as their cross- 
correlation. We also show computations of the covariances and Fisher matrices for these estimators, which will 
be needed in the subsequent statistical analysis. 

In order to refer to a specific case, in this work we adopted the guidelines set up for the planned Euclid 
weak lensing survey, as they are stated in the Euclid Red Book [4]. Different specifications are provided for 
different survey configurations. Here we chose the configuration covering a solid angle of Q = 15,000 deg 2 
over the extragalactic sky (f&y = 0.364). In this case the source redshift distribution reads 




(4.7) 




(4.8) 
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Figure 1. Left panel. The convergence power spectrum as a function of the angular scale for the fiducial Gaussian cosmo- 
logical model. The solid red line represents the total spectrum, while the two black dashed lines refer to the 1 -halo and 
2-halo contributions. Right panel. The convergence power spectra computed for the local shape non-Gaussian bispectrum 
with different values of the parameter /nl, as labeled. All spectra are divided by the convergence power spectrum in the 
reference Gaussian cosmology. 

with z* = 0.6374, a = 2, and B — 1.5. The corresponding median source redshift is z m = 0.9. 

As can be seen from the Fisher matrices reported in Eqs. (A. 51), (A.52), and (A. 53), a number of other 
parameters need to be specified in order to perform a statistical analysis. The first one is the average number 
density of background sources for which a shape measurement is feasible. The Euclid minimum requirement 
is « = 30 arcmin -2 , while its goal is « = 40 arcmin -2 . We adopted the former in this work. Also, we have 
to select the width of the multipole band where averages of the convergence power spectrum estimators are 
taken (see Eq. (A.39) and the discussion that follows). Consistently with the majority of previous works, we 
used A{ — 1 hereafter. The most delicate parameters to be selected are probably the intrinsic shape noises for 
the shear of background sources cr y , for the flexion cr^, and for their cross-correlation o~ x (see Appendix A for 
formal definitions), since they enter squared in the respective Fisher matrices. The value of the intrinsic shear 
noise suggested for Euclid is cr y = 0.3 (see however [94] for a different choice), but there are no standard 
values for cr^ and cr x . As a reference value we adopted the one suggested by [95], that is cr^ = 6.19 x 10 . 
This choice is discussed in more detail further below, where we also show the consequences of adopting a 
substantially smaller value. Lacking any evidence to the contrary, here and in the following we considered 
o~x = (CyO - !/)) 1 ^ 2 - This implies cr x = 43.1 for our current choice of parameters. With these assumptions, and 
considering the different scale dependencies of the flexion shot noise and the shear shot noise (compare Eq. 
A.40 with Eq. A.42), the former becomes smaller than the latter for I > 20, 000 [96]. 

5 Constraints on PNG 

5.1 Modifications to the convergence power spectrum 

In this Subsection we summarize the impact of PNG on the convergence power spectrum. Further details on this 
topic can be found in [37] . For reference, in the left panel of Figure 1 we show the convergence power spectrum 
evaluated in the reference Gaussian cosmology as a function of angular scale. We report the contributions of 
the 1-halo and the 2-halo terms of the non-linear matter power spectrum. As expected, the former dominates 
on very small angular scales, while the latter is more important at very large scales. Our subsequent analysis 
shall be limited at multipoles above 4nin = 50 (9 < 430 arcmin). While the Euclid Red Book recommends the 
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value £ m j„ — 5, we preferred to adopt the former value in order for the Limber approximation used above to be 
valid with high accuracy. 

In the right panel of Figure 1 we report the impact of PNG on the convergence power spectrum. For 
illustrative purposes we only consider non-Gaussianity with local shape, and a range of /nl values that brackets 
the current best constraints. As can be seen, the convergence power spectrum shows an increase at intermediate 
angular scales for positive skewness, and a decrease for negative ones. The peak of the difference is reached 
at ~ 30 arcmin angular scale, and it is somewhat larger than ~ 1% for /nl = ±100. Since the shear shot noise 
is white, while those of flexion and of their cross-correlation decrease toward small scales (see the covariances 
computed in Eqs. A. 40, A.42, and A.45), it is expected that flexion constraints on PNG would improve much 
faster than shear constraints when including smaller and smaller scales. In the next Subsection we show how 
this expectation is verified. 

Note that the behavior depicted in the right panel of Figure 1 is, for small angular scales, different from 
the results presented in, e.g., [97]. That work is based on numerical n-body simulations, and the modifications 
to the cosmic shear correlation function do not show any crossover. This is ascribable to the effect of PNG on 
the inner structure of dark matter halos, which we ignored, and makes our subsequent analysis conservative 
modulo the effect of baryons. 

5.2 Fisher matrix analysis 

We performed a Fisher matrix analysis in order to forecast the constraints that the convergence power spectrum, 
estimated either via shear, flexion, or their cross-correlation, can put on the level of PNG, depending on the 
shape of the primordial bispectrum. As already mentioned, the detailed calculations for these Fisher matrices 
are collected in Appendix A, to which we refer the interested reader. The fiducial cosmology that we adopted 
for this analysis is the standard ACDM model detailed in the Introduction, having /nl = 0. 

In this work we limited ourselves to study the level of non-Gaussianity /nl, and the least well-known of 
the standard cosmological parameters, erg. By the time Euclid will fly, it is reasonable to expect that all the 
other cosmological parameters will be known with great accuracy thanks to CMB experiments, such as Planck, 
in combination with Ho measurements. Hence it is fair to fix them to their fiducial values for our purposes. One 
note of caution with respect to the Fisher matrix analysis is in order. Here we assumed that the convergence 
power spectrum estimates obtained through shear, flexion, and their cross-correlation are independent of each 
other, so that the respective Fisher matrices can be simply summed. This assumption is acceptable for what 
concerns shear and flexion estimators, since they are measured from different moments of the galaxy light 
distributions. Still, when the cross-correlation between shear and flexion becomes relevant there might be 
other contributions to the Fisher matrix that have been neglected here [98]. Further studies will be needed on 
this point. 

In Figure 2 we show the joint constraints on the <x 8 - /nl plane for different CLs, given by the combi- 
nation of the shear-based estimate of the convergence power spectrum, the flexion-based estimate, and their 
cross-correlation. In Table 1 we report the lcr marginalized errors on the two individual parameters /nl and 
<x 8 obtained by using the same data. The marginalized lcr constraint for the j— th parameter can be found by 
inverting the Fisher matrix and then computing (Fjj) 1/2 - While £ m { a = 50 is left fixed, we adopted three differ- 
ent values for the maximum multipole £ max . In the first configuration we limited ourselves to relatively small 
multipoles, 6 max = 1,000, where the behavior of the matter power spectrum is sufficiently well understood, 
either thanks to numerical simulations or perturbation theory. In the second case we selected £ max = 3,000, 
which is about the upper limit within which the non-Gaussian part of the covariance (at least for the shear-based 
estimator of the convergence) can be neglected [99-102], as can the inner structure of dark matter halos and 
the impact of baryonic physics [91, 103, 104]. Finally, in the third situation we pushed the sum in Eqs. (A. 51), 
(A. 52), and (A.53) up to ^ max = 20, 000. This corresponds to an angular scale of 6 ~ 1 arcmin, and hence 
equals at considering the optimal situation in which the entire field covered by Euclid is tesselated in pixels of 
~ 1 arcmin 2 and each is usable for weak lensing studies. Although this latter configuration has been used by 
some authors (even recently, see [38]), it is quite extreme. A fruitful exploitation of such small scales would 
require detailed knowledge of the inner structure of dark matter halos, the impact of gas cooling, star formation, 
and AGN feedback, as well as the non-Gaussian part of the covariance. The physics of baryons, in particular, is 
expected to have a substantial impact at these small scales [104, 105], possibly erasing part of the non-Gaussian 
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Figure 2. Joint confidence regions for the parameters tr 8 and / NL (local bispectrum shape) given by the combination of 
cosmic shear, cosmic flexion, and their cross-correlation. The blue regions refer to 68.3% CL, the violet ones to 95.4% CL, 
and the pink ones to 99.7% CL. The primordial bispectrum shape is assumed to be local. The red stars at the intersection 
of the red dotted line show the fiducial Gaussian cosmology. In the different panels we show different combinations of the 
maximum multipole £ raax over which Fisher matrices are computed, and the intrinsic shape noise for flexion measurements 
cr„, as labeled. 



signal. As a consequence, the results for this configuration must be taken as rough approximations of the real 
situation. 

It should be noted that the multipole interval recommended by the Euclid Red Book ranges up to { max = 
5, 000. In our analysis however we adopted as a fiducial configuration £ lmx = 3, 000 for two reasons: first, we 
wanted to be sure that the non-Gaussian parts of the covariances are negligible, an assumption that loses more 
and more accuracy the smaller is the angular scale considered. Second, this allows a more direct comparison 
with our previous results in [37], where the exact same multipole interval has been used. 

Figure 2 and Table 1 show that by using the cosmic shear alone it is possible to constrain erg at the 
level of ~ 3.7 x 10~ 3 and /nl at the level of ~ 120 when £ max = 1,000. These constraints are only mildly 
improved when considering smaller scales, being reduced by ~ 20% if £ max = 3, 000 and by a factor of ~ 2 if 
Anax = 20, 000. The constraints on /nl are a factor of a few looser than those found in a previous work of ours, 
[37]. The reason for this is that at the time the specifications for Euclid were still not well defined, hence we 
adopted a background galaxy number density ~ 33% higher and an intrinsic shear noise ~ 36% smaller than 
in the present work. This resulted in a shot noise ~ 2.5 times smaller than here. Additionally, the assumed sky 
coverage was also one third larger there than here, and the adopted source redshift distributions are different. 

The constraints on cosmological parameters coming from flexion alone are much weaker than those due to 
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Table 1. The lcr marginalized errors on erg and /nl for the various estimators of the convergence power spectrum consid- 
ered in this work and their overall combination. The primordial non-Gaussian bispectrum is assumed to have local shape, 
and different combinations of i mia and have been considered (see the text for details). 
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Figure 3. The 68.3% CL joint constraints on erg and the level of PNG /nl (local shape) given by cosmic flexion (pink 
contours), cosmic shear (blue), and their cross-correlation (violet). The shape of primordial bispectrum is assumed to be 
local. In the different panels we show different combinations of the maximum multipole £ m!a over which Fisher matrices 
are computed, and the intrinsic shape noise for flexion measurements a v . Note that in the bottom right panel the contours 
stemming from flexion and the shear- flexion cross-correlation are not visible, being much smaller than the shear contour. 

shear alone. For £ max = 1 , 000 the forecasted error on cr% is of a few, and that on /nl is ~ 10 5 . Since we assumed 
that the number density of background sources usable for shape measurement be the same in the two cases, this 
difference is obviously due to the intrinsic shape noise, that for flexion is about 4 orders of magnitude larger 
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Table 2. The correlation coefficients between erg and /nl for the various convergence power spectrum estimators considered 
in this work and their full combination. A local shape for the PNG bispectrum is assumed. Different combinations of £ max 
and <t v are considered (see the text for details). 



o> = 6.19xl0 3 ov = 1.16xl0 2 

fmax = 1 . 000 f max = 3, 000 f max = 20, 000 f m ax = 3, 000 

cosmic shear -0.978 -0.987 -0.952 -0.987 

cosmic flexion -0.995 -0.992 0.492 -0.990 

cross-correlation -0.962 -0.985 -0.647 -0.990 

total -0.977 -0.987 -0.945 -0.990 



than for the shear. At the same time however, the cosmic shear shot noise is white, while the cosmic flexion 
shot noise decreases toward small angular scales. As a consequence, constraints on cosmological parameters 
from flexion alone are improved by a factor of ~ 6-7 for £ max = 3, 000 and by more than two orders of 
magnitude when C max — 20,000. This fact can be found in Figure 3, where we compare the 68.3% CL joint 
constraints on erg and /nl derived individually from cosmic shear, cosmic flexion, and their cross-correlation. 
Still, the constraints coming from cosmic flexion alone remain about one order of magnitude worse than those 
coming from shear alone. In order to have comparable constraints it would be necessary to extend the sum 
in the Fisher matrices substantially above the multipole for which the two statistical errors are comparable, 
£-20,000. 

The constraints on cosmological parameters given by the cross-correlation of shear and flexion are in 
between the two above. In particular, said constraints improve with decreasing scale faster than those due to 
shear alone but slower than those due to flexion alone. In the ideal configuration having £ max = 20, 000, the 
constraints are a factor of ~ 2 - 3 better than those coming from flexion, but still a factor of ~ 3 - 4 worse 
than those coming from shear. As a consequence of the above discussion, the constraints on cosmological 
parameters given by shear are basically unchanged by the inclusion of both the flexion and the cross-correlation 
between shear and flexion. The exception to this is given by the configuration with £ max = 20, 000, for which 
constraints are improved by ~ 8% upon inclusion of the flexion information. 

The improvement on cosmological constraints given by the inclusion of flexion for £ max = 20, 000 is 
not only due to the decrement of the statistical noise of the latter with increasing multipole, but also to another 
factor that can be better appreciated by examining Table 2. This table reports the correlation coefficients for the 
two cosmological parameters considered here, obtained by adopting shear, flexion, and their cross-correlation 
individually, as well as their full combination. In general, the correlation coefficient between the j— th and the 
k-th parameters is given by 

F~l 

r jk = . 1 ■ (5.1) 

As can be seen, in almost all circumstances the correlation coefficient between cr 8 and /nl is very close to -1, 
meaning that the two parameters at hand are almost perfectly anti-correlated with respect to the cosmological 
probes considered here. The correlation for cosmic flexion is slightly smaller than for cosmic shear, however 
the difference is only at the level of 1 -2%. This fact is of course expected because both the shear and the flexion 
covariances (as well as the covariance of their cross-correlation) depend in the same way on cosmological 
parameters. An exception is given by the configuration with £ max = 20, 000. In this case the correlation for 
cosmic shear remains quite close to -1, while that for cosmic flexion changes sign and becomes positive, 
taking a value of ~ 0.5. The correlation coefficient for the cross-correlation also shows a substantial increase, 
although it remains on the negative side. These changes can be better appreciated by looking at the bottom left 
panel of Figure 3. 

The reason for this behavior is easily understood by looking at Figure 1. At very small angular scales the 
effect of PNG on the convergence power spectrum gets reversed, with a positive /nl value implying a reduction 
of power and vice-versa. Since both the estimators of the convergence power spectrum based on flexion and on 
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Figure 4. Joint constraints on the level of non-Gaussianity / NL and cr 8 given by the combination of cosmic shear, cosmic 
flexion, and their cross-correlation. Pink contours refer to 68.3% CL, violet ones to 95.4% CL, and blue ones to 99.7% 
CL. Red stars at the intersections of red dotted lines mark the fiducial Gaussian cosmological model. Calculations assume 
''max = 3, 000 and a 9 = 6. 19 x 10 3 . Each panel refer to a different bispectrum shape, as labeled in the plot (see the text for 
more details). 

the cross-correlation between shear and flexion are weighted toward high multipoles, they are more strongly 
affected by this behavior with respect to the shear-based estimator. It follows that the improvement in the 
constraints given by the inclusion of flexion over those of shear alone are partially due to the degeneracy- 
breaking impact of flexion. We stress again that this result should be taken with caution. Even ignoring the non- 
Gaussian part of the covariance, the improvement due to the inclusion of flexion could be either overestimated 
or underestimated: the fact that dark matter halos are more concentrated in non-Gaussian cosmologies with 
positive skewness would decrease the scale at which the impact of /nl on the convergence power spectrum 
is reversed, hence reducing the improvement due to the inclusion of cosmic flexion. On the other hand, PNG 
might also imply a more intense formation of first stars [ 106, 107], with the subsequent larger energy feedback 
which would effectively blow up the inner parts of dark matter halos and/or avoid gas condensation at the rates 
of a ACDM universe, and hence work in the opposite direction. 

Let us now study in more detail the effect of the intrinsic flexion noise on the forecasting power of 
cosmological weak lensing. The value of cr v that we adopted in the preceding discussion was considered 
standard by [95]. At the same time, authors in [96] adopted a value that is slightly larger, while other authors 
[108, 109] considered numbers that are a factor of ~ 3 smaller than that considered here. Despite the fact that 
the values of cr^ used in all these work share the same order of magnitude, it is worth stressing that a standard 
value for the flexion intrinsic noise does not exist yet. The values selected in those of the aforementioned 
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works that attempt an actual measurement of cr^ both neglect the shear-flexion cross-talk [98] and make use 
of non-optimized weight functions for flexion measurements, hence they cannot be considered as thoroughly 
reliable. More specifically, the first flexion that we are considering here induces arc-like distortions on the 
images of background galaxies. Yet, galaxies are intrinsically very rarely arc-shaped, so a value of <x^ as high 
as the one adopted above might be unlikely. While pixel noise and measurement noise can indeed introduce 
arc-shaped distortions in galaxy images, these are systematic effects that are beyond the scope of the present 
statistical analysis. 

In order to estimate the impact of this uncertainty on cosmological constraints and to bracket viable 
alternatives, we repeated the Fisher matrix analysis for the configuration with ^ max = 3, 000, by replacing the 
cr^ used above with a substantially lower value. We set this value by requiring that the flexion S/N matches 
the shear S/N at the smallest scale in the configuration at hand, that is I = 3, 000. For a Singular Isothermal 
Sphere (SIS) the ratio between the flexion field and the shear field is 1/0, where is the angular separation 
from the center of the sphere (see, e.g., [110]). In the very outskirts of the SIS, say = 1 arcmin, this ratio 
equals 1/(60 arcsec). For a typical image 1 arcsec across, the ratio of the flexion signal to the shear signal is 
hence 1 /60 = 0.0167. Assuming that this holds also for LSS lensing and requiring that the flexion S/N and the 
shear S/N be the same at£- 3,000 implies cr^ = 387<x r 116. Although this choice admittedly bears some 
degree of arbitrariness, we feel it is the best that can be done until more detailed work on the topic at hand is 
performed. Moreover, this option also allows us to understand what level of cr^ would be necessary in order to 
obtain significant constraints on /nl from cosmic flexion. 

The results for the joint constraints on cr 8 and /nl implied by the new value of the intrinsic flexion noise 
are reported in the last two columns of Table 1, while the relative correlation coefficients are displayed in the 
last column of Table 2. As can be seen by inspecting these values, constraints due to cosmic flexion alone 
are improved by more than two orders of magnitude, while those due to the shear- flexion cross-correlation are 
improved by more than one order of magnitude. More importantly, both constraints are now tighter than those 
given by cosmic shear alone, by up to a factor of ~ 2 for cosmic flexion. This result might look surprising 
at first, since the intrinsic shape noise for flexion is still more than two orders of magnitude larger than that 
for shear. However, due to the scale dependence of the flexion shot noise, the latter becomes smaller than the 
shear shot noise already at £ > 400, and since the sum extends up to £ = 3, 000 overall the flexion becomes the 
dominant contribution. As a result, the total combination of shear, flexion, and their cross-correlation improves 
constraints on PNG by a factor of ~ 3 with respect to the shear alone. 

The dramatic change in constraining power due to a reduction in the flexion intrinsic noise can be appreci- 
ated also by looking at the bottom right panels of Figures 2 and 3. In the former we show the joint constraints on 
cr g and /nl given by the overall combination of the three estimators. The confidence ellipses are substantially 
smaller even than those resulting from the configuration with £ max = 20, 000 but with a high value of cr^. In 
the latter we compare the joint constraints given by each individual estimator taken separately. The difference 
in the sizes of the confidence ellipses with respect to the previously considered cases is so considerable that on 
the scale of the Figure only the shear constraints are visible, while the constraints stemming from flexion or 
the cross-correlation of shear and flexion are substantially smaller than those. Note the stark difference with 
the upper panels of the same Figure, where the confidence ellipse given by cosmic flexion alone did not even 
fit the scale. 

As a final step, we considered again the configuration with f max = 3,000 and cr^, = 6.19 x 10 3 and 
repeated the Fisher matrix analysis described above for three additional shapes of the PNG bispectrum, namely 
the equilateral, enfolded, and orthogonal shapes, as described in Section 2. The results for the joint constraints 
on crs and /nl obtained from the total combination of cosmic shear, cosmic flexion, and their cross-correlation 
are shown in Figure 4, while in Tables 3 and 4 we report the marginalized constraints for each individual probe 
and the correlation coefficients, respectively. As expected, forecasted errors on <r 8 are almost unchanged when 
considering different bispectrum shapes. On the other hand, constraints on the level of PNG are significantly 
looser for shapes different from the local one. This is so because the impact of non-Gaussianity on both the 
mass function and halo bias (and hence on the non-linear matter power spectrum) is weaker for these shapes 
(e.g., [ 14, 2 1 ]). Constraints on the level of PNG /nl range from ~ 340 for the equilateral shape up to ~ 500 for 
the orthogonal shape, a factor of up to five larger than the constraints for the local shape. 

Since we considered a high value of cr^, the inclusion of the flexion information does not bring any 
significant improvement on cosmological constraints for these non-Gaussian shapes, with the bounds on /nl 
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Table 3. The lcr marginalized errors on erg and /nl for the various convergence power spectrum estimators considered 
in this work and their overall combination. Four different bispectrum shapes are shown, as labeled. Calculations assume 
4n« = 3,000 and cr = 6.19 x 10 3 . 
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Table 4. The correlation coefficients between erg and /nl for the various convergence power spectrum estimators considered 
in this work and their overall combination. Four different bispectrum shapes for PNG are assumed, and calculations are 
performed with £ max = 3, 000 and <x v = 6.19 x 10 3 . 
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improving only by less than a percent upon combination of cosmic shear, cosmic flexion, and their cross- 
correlation. Also, the correlation coefficients are all very close to - 1 . The only exception is represented by 
the orthogonal model. In this model, the skewness of the density fluctuation distribution is well known to be 
positive for negative values of /nl and vice-versa, so that the direction of degeneracy between cr 8 and /nl 
is reversed with respect to other cases. This is very well visible in the bottom right panel of Figure 4. The 
forecasted constraints on the level of PNG obtained for bispectrum shapes other than the local one are roughly 
at the same level of the current bounds, derived either from the CMB or LSS tracers (see the discussion in 
Section 2). As explicitly shown for the local shape, a reduction of the flexion intrinsic shape noise would 
arguably improve these constraints by a substantial amount. In order to correctly address this issue however it 
is necessary to measure this intrinsic noise in a way that is both unbiased and optimized. Substantial work is 
thus still required in this direction. 

6 Discussion and conclusions 

In this work we analyzed the constraints on PNG that will be put by future wide-field weak lensing surveys such 
as Euclid. We considered shear, flexion, and their cross-correlation as estimators for the convergence power 
spectrum. In Appendix A we collected detailed calculations of the relevant covariances and Fisher matrices, 
which have been used in order to compute joint constraints for the level of PNG /nl and the amplitude of the 
matter power spectrum <x 8 using the accepted specifications for Euclid [4] . Our main results can be summarized 
as follows. 

• The latest Euclid specifications imply somewhat smaller sky coverage and average galaxy number den- 
sity than assumed in previous calculations. As a result, the constraining power of cosmic shear is weaker 
than formerly found. Specifically, the level of PNG for local shape bispectrum will be constrained at the 
level of A/nl ~ 100, which is at the same level of present-day constraints from the CMB and the LSS 
[63, 77]. 

• For the first time in the literature we forecasted joint Euclid constraints on /nl and <x 8 from cosmic 
shear for four different shapes of the primordial bispectrum, namely local, equilateral, enfolded, and 
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orthogonal. Constraints on the latter three non-Gaussian shapes are substantially looser than for the 
former, ranging between ~ 340 - 500, again at the level of current constraints from the CMB and the 
LSS. 

• We carefully reviewed the calculation of the covariance for a shear-based estimator of the convergence 
power spectrum, and used the same approach to compute, also for the first time, the covariances of a 
flexion-based estimator and of an estimator based on the cross-correlation of shear and flexion. These 
results have then been used to estimate the relevant Fisher matrices and hence combine the shear infor- 
mation with the flexion one. 

• The constraining power of cosmic flexion depends heavily on the assumed intrinsic shape noise <x^. 
Values adopted by several previous works might be too high to be realistic, and result in a constraining 
power for cosmic flexion alone that is much lower than for the cosmic shear: constraints on <x 8 are of 
order unity, while those on /nl are of the order of a few xlO 4 . The shear-flexion cross correlation is in 
between the two, with constraints still being barely significant with respect to cosmic shear alone. 

• Since the shot noise for cosmic shear is white, while that for cosmic flexion and for their cross-correlation 
decreases with decreasing angular scale, the higher the multipole considered in the Fisher matrix analy- 
sis, the more the constraining power of flexion grows with respect to that of the shear. Pushing the sum 
up to f max = 20, 000 improves the total joint constraints on /ml and erg by about ~ 8% with respect to the 
cosmic shear alone. Considering even larger multipoles would increase this relative contribution even 
more. 

• The value of <x^ suggested in the literature is likely neither unbiased nor optimal. A very rough order-of- 
magnitude estimate suggests a value for <r^ a factor of ~ 50 smaller than that, which produces an enor- 
mous improvement on the cosmological constraints due to flexion-related estimators. For local shape 
PNG flexion bounds become a factor of ~ 2 better than the shear bounds, and the overall combination of 
the three estimators performs a factor of ~ 3 better. 

As mentioned above, constraints on PNG obtained by setting ( max = 3, 000 are comparable with the 
current bounds derived from the CMB and the LSS [39, 63]. This is true for all the four bispectrum shapes that 
have been considered in this work. Only pushing the sum in the Fisher matrix analysis up to £ m - dx = 20, 000 
can produce forecasted errors that are substantially lower than the current constraints. However, we stress once 
more that these results should be taken with extreme caution. Firstly because we ignored the non-Gaussian part 
of the covariance in our analysis, stemming from the non-linear mixing of modes, which might be important 
at very small scales. Secondly because the non-linear clustering of matter and the effect of baryonic physics 
at small scales are still very uncertain (for a recent discussion see [111]). This consideration suggests several 
possible lines for future investigation. 

The non-Gaussian contribution to the cosmic shear covariance has been estimated in [102], and using 
a similar approach it should be possible to do the same for cosmic flexion and for their cross-correlation. 
Substantial effort is already being put in order to gauge the impact of gas cooling, star formation, and AGN 
feedback on the small-scale clustering of matter [91, 104, 105]. Although at the moment substantially scattered 
results are obtained due to the different implementations of non-gravitational physics, it is arguable that this 
issue will be much better under control by the time Euclid data will be available. Last but not least, it is 
important to better understand the impact of PNG on the internal structure of dark matter halos. While it 
is relatively certain that halos in non-Gaussian cosmologies with positive skewness will be more compact 
than their Gaussian counterparts, there still are many obscure points. Just to mention one, nobody has yet 
investigated if indeed the NFW shape is such a good fit to the density profiles of dark matter halos in models 
with PNG as it is for Gaussian cosmologies. 

The results of this paper show that if the flexion intrinsic shape noise is at the level quoted by several 
previous works, cosmic flexion alone will not be usable for constraining cosmological parameters, not even 
with future state-of-the-art space missions, unless a better exploitment of highly non-linear scales becomes 
feasible. If cr v turns out to be smaller instead, constraints due to cosmic flexion can be comparable or better 
than those due to cosmic shear. In the latter case, flexion will be a very important addition to the weak lensing 
analysis, and will allow to substantially improve constraints on PNG. In the former, the contribution of cosmic 
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flexion to the constraining power of cosmological weak lensing can be safely neglected. Either way, it should 
be borne in mind that the final Euclid data product will not be limited to weak lensing alone, and will overall 
allow to constrain the level of PNG at the level of A/nl ~ a few [78]. 
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A Covariances and Fisher Matrices for Cosmological Weak Lensing 

In this Appendix we describe in detail how we computed the covariances of cosmic shear, cosmic flexion, and 
their cross-correlation, as well as the related Fisher matrices. First, let us recall how the shear is related to the 
convergence in Fourier space, namely 



which motivates us to redefine the shear in Fourier space simply as y(t) — k{{). In other words, we are adopting 
a scalar representation for the shear tensor field. 

Similar considerations can be applied to the flexion. We shall limit ourselves to the first flexion <p(0), 
leaving the second flexion for future studies. First, let us recall the definition of the complex flexion, <p(Q) = 
8k{6) — 8\k{6) + \82k{6). Then, from the properties of the Fourier transform it easily follows that 




(A.l) 



From this relation it easily follows that 



(A.2) 



(A3) 



so that 
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= k{t)k*{e) (e\ + 4) = e 2 k(w*(t) ■ (A.4) 

Hence we redefine the flexion in Fourier space as (p({) = I k{t) (we might call this the absolute flexion, as 
it loses the directional information). We recall that the convergence is a real-valued field, hence its Fourier 
transform enjoys the symmetry k*{{) = k(-t). The redefined shear and flexion also share the same symmetry, 
as it is easy to verify. 

By using the Limber approximation the power spectrum of the convergence can be written as 

(k(f)k*(n) = {2n) 2 5 D {t - t')P K {t) , (A.5) 

where the angular brackets represent the ensemble average and P K {t) = t^P^it), with P<f,(f) being the power 
spectrum of the lensing potential <p (see also Eq. 4.7). By using the previous relations, the power spectrum of 
the flexion can now be written as 

{(pdWiO) = (2n) 2 S D (e - t')P^{{) = ££'(2K) 2 6 D (e - e')P K (0 ■ (A.6) 

It follows that the flexion power spectrum can be related to the convergence power spectrum by P^(f) = ( 2 P K (t). 
This is a well known relation obtained for the first time by [9] . Analogously, the cross-power spectrum between 
the convergence and the flexion reads 

| mWiO) + {(p{t)K it'))] = (2n) 2 6 D (e - e')P x (0 = \tf + O (2n) 2 5 D {t - f)P K {f) , (A.7) 

which implies that the cross correlation between the convergence and the flexion can be related to the conver- 
gence power spectrum as P x {t) = i P K (f). This also coincides with the result shown by [9]. 

A.l Estimators 

In order to compute the covariances, we need first to quantify estimators for the convergence power spectrum 
as measured through shear and through flexion. For this we shall follow closely the recent work by [96]. Let us 
suppose to have shape measurements (either ellipticity for the shear or some other moment of the luminosity 
distribution for the flexion) for a sample of N galaxies uniformly distributed across a field of area A-nf^y, so that 
the average number density of galaxies is h = NI{Ajif^ y ). According to the authors in [112, 113], an estimator 
for the complex shear in Fourier space can be written as y (e \l) = h y{t) + s y ({), where s y (8) represents 
the intrinsic contribution to the measured ellipticity, and since this can be measured only at an actual galaxy 
position, its Fourier transform can be discretized according to 

N 

s y (f) = ^ s 7 {0j) exp {-it ■ 0j) . (A.8) 

7=1 

We further recall that the scalar representation of the shear reads y{{) = k(t), hence an estimator for the 
convergence based on shear measurements can be written as 

k y {t) = f (e \f) = n k(f) + s 7 (f> . (A.9) 

In a completely analogous way, we can come up with an estimator for the flexion, £> <e) (f) = n (p{t) + e^(f), 
where the function e^ff) encapsulates this time the intrinsic contribution to the measured flexion. Furthermore, 
because of the way we redefined the flexion here, we know that (p(() = I k(t), hence an estimator for the 
convergence based on flexion measurements can be written as 

ky(f)=nm + -e lp {t). (A. 10) 

Note that, since the intrinsic shear and flexion are assumed to be randomly distributed, then (s y (0j)) — = 
(e v (0j)) for all 

We now introduce a shear-based estimator for the convergence power spectrum as [113] 
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(A. 11) 



47r/ sky « 2 j nf 

where Q.( is some small region in multipole space around I, m(D.() is its volume, and f^ y is the fraction of sky 
covered by the shear (flexion later on) measurement [114]. The quantity <x y represents the average dispersion 
of the intrinsic source ellipticity, that is defined as 



(s y (ej)s;(e k )) = s*(r. 

The ensemble average of this estimator reads 

(2 



(A. 12) 



= T-h-2 f 4^ <WW) - 5 > (A.13) 
where the argument of the integral in the previous equation can be estimated thanks to the following reasoning, 



(k y {t)k;{t')) = h 2 (2n) 2 6 D (( - op At) + h (s r (e)k\e')) + n (ms;(0) + (s y {t)s y {t')) . (A.14) 

The second and the third terms of the previous sum represent the cross-correlation between the convergence 
field and the intrinsic ellipticity of sources. The two can be considered to be statistically independent, so that 
these two terms vanish. The last term is instead the power spectrum of intrinsic ellipticity, and it is a little bit 
trickier to determine. Specifically, 



{s y (t)s y (t')) 



J] s y (0j) exp (-it ■ 9j) J] s* 7 (0 k ) exp (if ■ k ) 

;=1 JU=1 



= 2 (sy(Oj)e;(0 k )) exp (-it ■ Oj + if ■ e k ) = 



= o* J] exp [-i{l - {') ■ 0j] = cr 2 y h (2n) 2 5 D {{ - f) 

7=1 



(A.15) 



By replacing this result in the previous Eq. (A.14) we obtain for the ensemble average of the shear-based 
estimator of the convergence power spectrum 



PM) + 



(A. 16) 



By assuming that the region Cl( is small enough, so that the function in square brackets does not vary signif- 
icantly across it, we can pull the function itself outside the integral. Moreover, the quantity (27t) 2 5d(0 can 
be interpreted as the Fourier transform of the window function for the survey at hand, as long as the average 
separation between neighboring galaxies is much smaller than the extent of the field [113]. When computed at 
zero, this just equals the area of the survey, 47r/ s ] £ y. It follows that 



(p ( J\e)) - p K {t) . 



(A.17) 



Hence the estimator for the convergence power spectrum defined in Eq. (A. 1 1) is unbiased. 

In a similar way, it is possible to define a flexion-based estimator for the convergence power spectrum, as 



4;r/ sky n 2 J n m{Q.{) Y * nt 1 



/sky" 

where in this case the average dispersion of the intrinsic flexion is defined by 



(A.18) 



-24- 



(e v (Oj)s;{e k )) = 5*0% 



(A. 19) 



Note that in this case the last term on the right-hand side of Eq. (A. 18) depends on the multipole, consequence 
of the fact that the noise for a flexion-based convergence estimate is not white [96]. Specifically, by considering 
the ensemble average of the estimator in Eq. (A. 18), we obtain 



(pf (/)) = f ^n(2*) 2 <5 D (0) 

x ' 47r/ sky J ai m(Llf) 



PM) + 



and with the same approximation used above about fif it follows that 



(/><%)) - P K (f> . 



(A.20) 



(A.21) 



Similar considerations can be applied in order to construct an estimator of the convergence power spec- 
trum that is based on the cross-correlation of shear and flexion. Such an estimator can be written as 



2 47r/ sky « 2 J n , m(Llt) L v r J 



where now 



(A.22) 



(A.23) 



The ensemble average of this estimator reads 

< P ^> 4^ 1 ^ [<WW> + <^«W>] - Yi • (A - 24) 

and adopting the same procedure highlighted above the argument of the integral in Eq. (A.22) can be explicitly 
recast by using the following expression, 



(k y (t)k;{t')) + {^m;(e')) = 2(2n) 2 6 D (e - n 



9 (rill \ 

h 2 P K {{) + «-i - + - 



2 \€ 



The last term of the previous equation simplifies when t = (', so that 



x 1 47r/ sky J nf m(D f ) 



PM) + 



ilC 



PM 



(A.25) 



(A.26) 



Once again, the convergence power spectrum estimate resulting from the cross-correlation of shear and flexion 
has a colored noise, which decreases with increasing multipole, albeit more slowly than the noise on the purely 
flexion-based estimate. 

A.2 Covariances 

From the discussion above, it is quite clear that the covariance of the convergence power spectrum will depend 
on the observable that is adopted for estimating it. Let us start with the covariance of the shear-based estimate, 



v [p ( i\e\ pf{e')\ = {p { At)P ( ? (O) - (H y) w) (H r V)) - {pf{€)P ( I\r)) - p K {t)PM) . (A.2?) 

The first term of the sum above requires knowledge of the convergence four point correlation function. Specif- 
ically, 
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x ' (47r/ sky ) 2 n 4 J nf m(Q^) J^, m(Q r ) \ 7 7 ' 



" FT^h3 f "777^ (WW)~ 
" . 7 -3 ~7FT^\W)W)/ + 



(A.28) 



According to the calculations performed in the previous Subsection, the integrals in the second and third terms 

(y) 

on the right-hand side of the previous equation can be computed similarly to the ensemble average of P K {I) 
and Pj\('), respectively, thus giving 

(p ( /\e)p7\n) = / , 2 _ 4 f -^fr f ^(WWW')W'))- 



+ 4 

n 



(A.29) 



It is well established [101, 115] that thanks to the Wick's theorem [116] the four point correlation function 
can be written as a combination of two point correlation functions, with an additional connected part remaining. 
In the specific case, 



(V0WW)W)) = (WW)(W)W>) + (WW) (WW)) + 

+ (*y(£)W>) (WW)) + (w W W) W)) c ■ (A30) 

The double integral of the first term of the sum above gives 

(47r/ sky ) 2 n 4 J nf m(D f ) Ja t , m{Q. e ) \ 7 ' \ 7 ' n n « 2 

(A.31) 

which cancel, respectively, with the second term of the covariance (last term in Eq. A. 27) and with the last 
three terms of Eq. (A.29). There remain three terms in the covariance, namely 

y[^W,^V)] = * 4 f 4n f 4fr(WW))(W)W')) + 

L J (4^/sky) 2 " 4 Jn, m(€lt) J Q[ , m(Q r ) x ' v 7 7 ' 

+ ^ / , 2 -4 f -777T f 4?f^(WW'))(WW')) + 
(4^/sky) 2 " 4 Jn f wCQf) J n ,, m(Q r ) x ' v ' 

+ TTT-^l f -77^ f 4^(WW)W)W')) C ■ (A.32) 
(47r/sky) 2 « 4 Jn f m(€l() Ja,, m{D.f) \ 7 y 'c 

The last term in the above sum is just an integral over the trispectmm of the estimated convergence, and arises 
due to the mode-coupling induced by the non-linear gravitational growth of structures. If one limits himself 
or herself to relatively large scales, the impact of this term can be neglected. We adopted this approximation 
in the present work, hence we simply label this integral as t 7 ((, {') and disregard it later on. In order to make 
progress, we recall that k{{) - k*(—f). Moreover we assume that the real and imaginary parts of s 7 (t) are 
statistically independent [113]. It easily follows that 

{s y {t)s y {t')) = (s 7 {t)e 7 {-t')) (A.33) 
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(despite being s y {0) still a complex-valued function), and eventually 



(k y (t)k y (n) = (ty(Wy(-0) ■ (A.34) 

Making use of the former assumption we can recast the covariance of shear-based estimates as 



L J (47r/ sky ) 2 Jn f m{Q. e ) J Qe , m(Ll e ) 

1 £ (~* 

rr-^ -t?^ -77T(2^ D (^-r)ci r) (^r)(2^D^-r)ci r) (r^) 



+ 



(4* 

+ Ty({, (') 



where 



C ( j\lJ') = P K {t)+^ • 



+ 

(A.35) 
(A.36) 



Note that the function C^\t , I') is actually independent of the second argument, however we left it explicitly 
indicated for consistency with the flexion-based estimate discussed below. Now, since the convergence power 
spectrum depends only on the modulus of the multipole vector, it follows that C ( J\l, (') = \—t, ('). Also, 
we assume that the integration area Of is left unchanged under the transformation t — > -t, so that the previous 
equation can be recast as 

^\p { j\t),P ( j\t')U 2 f f -^(27r) 2 fe^-r)Ci y) (^r)(2^D(^')C< r) (r^)+T y (^r) 

L J (47r/ sky )- J a, m{Q. e ) J a[l m(Qf) 

(A.37) 

It is easy to see that the above double integral does not vanish if and only if the two regions in multipole 
space Q.( and Clf coincide. Since we eventually considered circularly symmetric regions in multipole space 
and the modulus t of the multipole vector actually takes integer values, it follows that 



L J (4^/ s ky)- Ja, m z (Q. e ) J 



= 6 



By extracting the function out of the integral, we obtain 



(A.38) 



/ sky m(Qf) 



PM + -4- 

n 



+ T y (e,e'). (A.39) 



In many circumstances it is assumed that Q.( is just the shell in multipole space included between C and I + At, 
where it is required that Alji <sc 1. If this is the case, then m(Q f ) = 2jt£A£ + nAt 2 = nAt(2£ + AC)- It follows 
that the covariance can be written as 



At(2( + A^)/ Sky 



P K (0 + -=■ 



(A.40) 



This corresponds, modulo the non-Gaussian part, to the error expression first derived by [92] and used, e.g., by 
[117]. 

A similar procedure and the same approximations can be used for the flexion-based estimate of the con- 
vergence power spectrum. The only relevant difference consists in replacing the function C^\t, I') with 
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cf\{J') = P K {0 + 



so that, eventually, 



where 



M(2( + AO/sky 



hW 



PM + 



(A.41) 



(A.42) 



(A.43) 



is the contribution to the covariance coming from the non-linear clustering of matter. 

It now remains to be calculated the covariance of the convergence power spectrum estimator based on the 
cross-correlation of shear and flexion. The details in this case are more cumbersome to work out, hence we 
report here only the most important steps and the final result. As usual, the covariance can be defined as 



-§ PM ~^ p ^-^- PMP ^- (A - 44) 



The ensemble average inside the double integral above returns the sum of four ensemble averages, that thanks 
to the Wick's theorem provide a total of sixteen terms. The terms that do not mix functions depending on £ 
with those depending on ij' can be rearranged to cancel the last row of Eq. ( A.44). The remaining terms can be 
collected together by using the assumptions above, and return the following final expression: 



Pk{D + -=- 
n 



p k (0 + ^ 

nl 



+ T x {tJ'), (A.45) 



where 



+ (k v (^k;(^k 7 ^')k;(f)) c + (^(OWVf ))J (A.46) 

encapsulates the contribution due to the non-linear clustering of matter. 
A.3 Fisher matrices 

In order to compute the Fisher matrix for a generic estimator of the convergence power spectrum, let us assume 
that such estimator is measured for all multipole values included between 6 m i n and £ max , and let us leave A£ 
unspecified. Under the assumption of a multivariate Gaussian likelihood the Fisher matrix can then be derived 
from the covariance as 

Fjk = ^Tr [Aj(t, f)A k {t, O + (/ , e'Wjktf, (')] , (A.47) 
where the moduli of both t and I ' run from to i max , and 

AM, t') = ^{.t, f )— *¥(€, {') , (A.48) 
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while 



A#;t(/, O = ^-PAD^-PAO + ^-PAD^-PAO ■ (A.49) 

The vector ^ represents the set of cosmological parameters on which the convergence power spectrum is as- 
sumed to depend on. 

In order to be more specific, let us write down the Fisher matrices for the various convergence estimates 
considered in this Appendix. The Fisher matrix for the shear-based estimate easily follows from 

M , 26 k d 

Af{t, t') = — a-PM • (A.50) 

1 P K {t) + rf/n dqj 

It should be noted that the trace of the matrix in Eq. (A. 47) is just a sum over the diagonal terms, for which 
t — {'. It follows that the Fisher matrix can be simplified to 

F 1k= Y, — r — -ir-ir p M)Tr p M ■ (A - 51) 

[PAD + rf/n] d 1J ^ 

A similar expression holds for the Fisher matrix related to the flexion-based estimate of the convergence power 
spectrum, namely 



F 



7 " -^^PAt)—P K {t)- (A.52) 

An expression that, for high multipoles, is equivalent to this, has been derived along a different route by [118]. 

The Fisher matrix for the convergence power spectrum estimator based on the cross-correlation of shear 
and flexion is a little bit more cumbersome, but straightforward to compute nevertheless. The final result can 
be written as 



= V 7 tt ; \2A{(2£ + M)f sky 

e tii [PAD + o*/n] [PAD + crl/(n{ 2 )] + [PAD + o*l(ffi? [ 

16[PAD + (r 2 7 /(4n) + crl/(4h( 2 ) + C r 2 x /(2hD] 2 \ d 3 



—PADt-PAD ■ (A.53) 



[pad + vj/n] [pad + (rime 2 )] + [pad + rt/m] 2 J d( U dc & 

These Fisher matrices can now be combined together to provide forecasted constraints on PNG. 
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